function fv=ML2(y1,x,r0,a,b,runmax,dmy,tolv,tau,y2r0,yir0,frac)
%******************************************************
%** This program is to execute maximum likelihood    **
%** estimation for ML2xtreg.m                        **
%** Created by Min-Hua Huang, Jan 25, 2012           **
%** Model: Maximum Likelihhod by Fixing at Minimum,  **
%**        with Constrained Optimization, Conducting ** 
%**        Correction on Location Estimates.         **
%** y1: inividual-level dependent variable           **
%** x: covariate matrix                              **
%** r0: initial parameter value                      **
%** a: lower bound                                   **
%** b: upper bound                                   **
%** runmax: number of maximum iterations             **
%** dmy: number of dummy variables                   **
%** tolv: tolerance value                            **
%** tau: step-size factor                            **
%** y2r0: grandmu and sd                             **
%** yir0: groupmu and sd                             **
%** frac: rescaled parameter, eta                    **
%******************************************************
format long; warning off all;
% basic setups
m=size(x,2)-1; record=ones(1,4*m+17)*-999; 
hr=boundary1(y2r0,yir0,frac); b=b-hr(1); a=a-hr(2);
% initial values
lamda=ones(2*m+6,1)*10; kk=1; run=1; stop=1; n=size(y1,1);
%****************** Record the initial results ***************************
sumD1=0;
for i=1:1:n
        if abs(b-x(i,:)*r0(1:m+1))<5*r0(m+2) || abs(a-x(i,:)*r0(1:m+1))<5*r0(m+2)
           D=normcdf((b-x(i,:)*r0(1:m+1))/r0(m+2))-normcdf((a-x(i,:)*r0(1:m+1))/r0(m+2));           
        else
           D=1; 
        end   
  sumD1=sumD1-log(D*sqrt(2*pi)*r0(m+2));
end
record(1,1)=0;
record(1,2:m+3)=r0';
record(1,m+4:2*m+5)=zeros(1,m+2);
record(1,2*m+6:4*m+11)=lamda';
record(1,4*m+12)=sumD1-1/(2*r0(m+2)^2)*((y1-x*r0(1:m+1))'*(y1-x*r0(1:m+1)));
record(1,4*m+13:4*m+17)=ones(1,5)*-999;
% maximum and miminum vector of covariates 
xmax=zeros(m+1,1); xmin=zeros(m+1,1);
   for i=1:1:m+1
       xmax(i)=max(x(:,i));
       xmin(i)=min(x(:,i));
   end    
xlimit=[xmax xmin];   
fprintf('Iteration % 4.0f\n',kk);  
while stop>0 && run<runmax
% vectors of indicator variables
v=vsign(r0);
% constraint vector
ci=cpar(r0,a,b,x,xlimit,v,dmy);
% gradient vector
%df=dfpar(r0,a,b,y,x); 
%L=Lpar(r0,a,b,y,x,tau);
overall=DH(r0,a,b,y1,x,tau);
df=overall(1,:)';
L=overall(2:m+3,1:m+2);
% first derivative of the constraint vector
ai=aipar(r0,x,xlimit,v,dmy); 
%****************** Step 1:check KKT conditions ********************
% KKT
KKT1=df+ai'*lamda; check1=sum(abs(KKT1)>tolv);
KKT2=ci; check2=sum(KKT2>tolv);
KKT3=lamda; check3=sum(KKT3<-tolv);
KKT4=lamda'*ci; check4=sum(abs(KKT4)>tolv);
check=check1+check2; %+check3+check4;
   if check==0
      stop=0; 
   end
% If KKT conditions are not satisfied, then iteration continues.
%****************** Step 2:compute dk and lamdaQP ********************
% Hessian matrix
% solving the QP problem  
opts = optimset('Display','off');
d=quadprog(L,df,ai,-ci,[],[],[],[],[],opts);  
lamdaQP=(ai')\(-df-L*d);
%****************** Step 3:set new x and lamda ***********************
r1=r0+d;
lamda=lamdaQP;
%****************** Record the updated results ***********************
sumD2=0;
  for i=1:1:n
        if abs(b-x(i,:)*r1(1:m+1))<5*r1(m+2) || abs(a-x(i,:)*r1(1:m+1))<5*r1(m+2)
           D=normcdf((b-x(i,:)*r1(1:m+1))/r1(m+2))-normcdf((a-x(i,:)*r1(1:m+1))/r1(m+2));           
        else
           D=1; 
        end   
  sumD2=sumD2-log(D*sqrt(2*pi)*r1(m+2));
  end  
record(kk+1,1)=kk;
record(kk+1,2:m+3)=r1';
record(kk+1,m+4:2*m+5)=d';
record(kk+1,2*m+6:4*m+11)=lamdaQP';
record(kk+1,4*m+12)=sumD2-1/(2*r1(m+2)^2)*(y1-x*r1(1:m+1))'*(y1-x*r1(1:m+1));
record(kk+1,4*m+13)=record(kk+1,4*m+12)-record(kk,4*m+12);
record(kk+1,4*m+14:4*m+17)=[check1 check2 check3 check4];
%****************** Step 4 and 5: return to Step 1 *******************
r0=r1;
kk=kk+1;
fprintf('Iteration % 4.0f\n',kk);    
run=run+1;
end
maxL=-10^20; maxt=-999;
for i=1:1:kk-1
    if record(i,4*m+12)>maxL && record(i+1,4*m+15)==0
       maxL=record(i,4*m+12);
       maxt=i; 
    end
end
if maxt~=-999 
   fv=[record(maxt,2:m+3) record(maxt+1,4*m+15) record(maxt,4*m+12)];
else    
   fv=ones(1,m+4)*-999; 
end   
  